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ABSTRACT 

We propose a new Monte Carlo method to study extended X-ray sources with the European Photon 
Imaging Camera (EPIC) aboard XMM-Newton. The smoothed particle inference (SPI) technique, 
described in a companion paper, is applied here to the EPIC data for the clusters of galaxies Abell 1689, 
Centaurus, and RX J0658-55 (the "bullet cluster"). We aim to show the advantages of this method 
of simultaneous spectral and spatial modeling over traditional X-ray spectral analysis. In Abell 1689 
we confirm our earlier findings about structure in the temperature distribution and produce a high- 
resolution temperature map. We also find a hint of velocity structure within the gas, consistent with 
previous findings. In the bullet cluster, RX J0658-55, we produce the highest resolution temperature 
map ever to be published of this cluster, allowing us to trace what looks like the trail of the motion 
of the bullet in the cluster. We even detect a south-to- north temperature gradient within the bullet 
itself. In the Centaurus cluster we detect, by dividing up the luminosity of the cluster in bands of 
gas temperatures, a striking feature to the north-east of the cluster core. We hypothesize that this 
feature is caused by a subcluster left over from a substantial merger that slightly displaced the core. 
We conclude that our method is very powerful in determining the spatial distributions of plasma 
temperatures and very useful for systematic studies in cluster structure. 

Subject headings: galaxies: clusters: individual (Abell 1689, IE 0657-55, Centaurus) — methods: data 
analysis — x-rays: galaxies: clusters 



1. INTRODUCTION 

Early work on X-ray emission from clusters of galax- 
ies, based on X-ray data obtained with instruments of 
modest angular and spectral resolution, implied that the 
profiles of X-ray emission are smooth and that the spec- 
tra can be adequately described as nearly isothermal hot 
plasma, generally indicating relaxed structure. This pic- 
ture has changed markedly with precise imaging data 
from the Chandra and XMM-Newton instruments: X- 
ray images of clusters reveal complex intensity distri- 
butions, where in general, the surface brightness lacks 
circular symmetry. The cluster emission cannot be de- 
scribed by a single- temperature plasma. Neither is the 
distribution of temperatures spherically symmetric, and 
it cannot simply be des cribed as radially dependent (e.g. 
iMarkevitch et al.ll2000l ). In addition, simple but well- 
motivated models such as "cooling flow s" fail to ade- 
quate ly describe the observations (e.g. iPeterson et al.l 
l200lh . even for clusters that are otherwise "relaxed." 
This is likely due to a complex history and physical pro- 
cesses associated with the cluster formation, which is yet 
to be fully understood. Such complexity may include 
effects of recent merger activity; large-scale "bubbles," 
presumably due to the interaction of the outflows pro- 
duced by the central active galaxy; sharp abundance gra- 
dients associated with recently trig gered star formation ; 
or other, still unknown processes (Sa nders et al.l (20041 ). 
Clearly, an analysis technique not relying on symmetry 
of the flux or temperature distribution is needed. 
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The methods considered for this task include fit- 
ting isothermal spectra to fixed or adaptively binned 
grids of photon s acr o ss the detect o r pla ne (e.g. 
IMarkevitch et all 120001 : ISanders et all l2004f ). An- 
other approach is imagi ng deprojecti on or "onion- 
peeling" methods (e.g. iFabian et al.1 Il98ll ). which 
in turn h ave been extended to spectroscop i c de- 
projection (lArabadiis et al.l 120021: lArnaud et al. 1200 ll : 



iKaastra et al.1 12QQ4J : lAndersson fc Madeiskil l20ol . An 



alternative method developed recently by some of 
us and termed "smoothed particl e infe rence" or 
SPI ([Peterson. Marshall, fc Anderssonl (|2007f ): hereafter 
PMA07) relies on a description of a cluster as a large set 
of "primitives," which in our case are smoothed parti- 
cles, or "blobs." Each of those is described by its overall 
luminosity and a spatial position, but also by a Gaus- 
sian width, a single temperature, and a set of elemental 
abundances. A large set (well upward of a hundred) of 
such primitives is then propagated through the instru- 
ment response using Monte Carlo (MC) techniques, and 
their parameters are adjusted via the use of Markov chain 
methods to map out the likelihood of the distribution as 
compared to the observation. The resulting distribution 
is then a good "nonparametric" description of the clus- 
ter. 

In this paper, we describe the implementation of this 
method to observations obtained with the XMM-Newton 
imaging detectors known collectively as the European 
Photon Imaging Camera (EPIC). We apply the method 
to data from three clusters of galaxies, namely Abell 
1689, RX J0658-55, and the Centaurus cluster. All 
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three clusters exhibit a substantial degree of complex- 
ity, and furthermore, all three have good-quality, deep 
XMM-Newton observations, which in turn can be com- 
pared against previous X-ray analyses as well as data 
from other instruments, such as the Chandra observa- 
tory. In particular, Abell 1689 is likely a merger or a 
superposition of two components aligned clo se to the line 
of sight (e.g. lAndersson fc Madeiskill2004ft : RX J0658- 
55 r eveals an ongoing merg er close to the plane of the 
sky ([Markevitch et aLll2002l ): and the Centaurus cluster 
shows abundance gradients as well as filaments and bub- 
bles, presumably caused by the ene rgy deposited in th e 
cluster by the central radio source ([Fabian et al.ll2005l ). 
The choice of objects should provide a good illustration 
(or test) of the method for three quite different cases. 
To perform this analysis, we constructed a Monte Carlo 
for the XMM-Newton EPIC detectors that is analogous 
to the MC that some of us developed for the Re flecting 
Grating Spectrometer (RGS; [Peterson et al.ll2QQ3 ). How- 
ever, here we use the new Markov chain-based method, 
SPI, for parameter iteration, as described above. 

Below, in Section 2, we describe the previous and cur- 
rent observations of the three clusters and summarize 
the data reduction procedures. In Section 3, we sum- 
marize the SPI method and describe the specifics of the 
XMM-Newton EPIC response function used here. We 
also describe the application of the method to the data. 
In Section 4 we discuss the results of the analysis, includ- 
ing the spatially resolved maps of the spectral parame- 
ters. In Section 5 we summarize the paper and discuss 
the advantages and disadvantages of the method and, fi- 
nally, we discuss the future applications of the method 
in Section 6. 

2. CHOICE OF TARGETS AND DATA REDUCTION 

Here we describe the choice of targets selected to 
demonstrate the capabilities and versatility of our 
method. The three chosen targets span a broad range of 
complexity of the flux, temperature, and redshift. Note 
that the exact same analysis chain is applied to all ob- 
jects. 

2.1. Abell 1689 

This cluster is one of the objects most studied with 
gravitational lensing techniques, as the optical im- 
ages clearly rev eal arcs and arclets, a llowing strong- 
lensing analysis ([Tyson fc Fischer! 1 1995f ). Deep studies 
with the ESO MPG Wide Field Imager provide addi- 
tional constraints towards the determination of the clus- 
ter gravitational potential, v i a the use of we ak lens- 
ing ([Clowe fc Schneider! I2QQ1I : IKing et all |2QQ2[ ) . More 
recently, the mass profile has been detailed further 
with combined strong and weak lensing, using Hub- 
ble Space Telescope (H ST) ACS and Subaru images 
([Broadhurst et all l2QQ5h . The optical data (including 
studies of galaxy veloci ties) indicate that the cluster 
contains sub-struc tures ([Miralda-Escude fc Babul|[l995l ; 
lLokas et al. 1120061 ) . The early X-ray data for this cluster 
indicated an X-ray intensity profile that implied a single, 
relaxed system, but the mass determination from the X- 
ray analysis indicated a lower mass (by about a factor 
of 2) than the mass value dete rmined from lensing (e.g. 
IMiralda-Escude Babullfl995h . 

The analysis of the XMM-Newton observation used 



in this paper, but using more traditional tech- 
niques than the SP I met hod, was presented by 
lAndersson fc Madeiskil (|2004 ) , and we refer the reader to 
that paper for more extensive discussion of Abell 1689's 
properties as well as previous X-ray observations. In 
summary, the X-ray emission appears to be symmet- 
ric, and the average temperature inferred for the re- 
gion of 3' in radius is 9.3 ± 0.2 keV, assuming the value 
of the Galactic column of 1.8 x 10 20 cm -2 derived by 
iDickev fc Lockmanl ([I990f ). However, the spatial analy- 
sis indicates that the temperature of the emitting plasma 
is clearly not uniform, ranging from ~ 7 to ~ 10 keV, 
with a hint of a temperature gradient in the southwest 
- northeast direction. In addition, the redshift of the 
emitting gas as inferred from the X-ray spectrum varies 
across the image, with a high-redshift structure to the 
east, with z = 0.185 ± 0.006, separated from the rest of 
the cluster at z = 0.17. This strongly suggests that the 
cluster consists of two components in projection, which 
either have started to merge, or are falling towards each 
other. One of the premises of our study is to determine 
if the two sub-components are indeed related, or if they 
should be treated as two separate clusters that happen 
to be located close to the same line of sight in the sky 
and are observed in projection against each other. 

2.2. RX J0658-5557 

RX J0658-5557 (or 1ES 0657-55.8), at z = 0.296 
was first disco vered as an extended X-ray source by the 
Einstein IPC ([Tucker et al. I fl995^> and was later found 
from Advance d Satellite for Cosm ology and Astrophysics 
(ASCA ) data ([Tucker et al .1119981 ) to have a temperature 
of about 17 keV, although subsequent simult aneous anal- 
ysis o f ASCA and ROSAT PSPC data by iLiang et all 
([2000D suggested a lower value of 14.5t 2 ;£ keV. Still, even 
the revised value of the temperature makes it one of the 
hottest known clusters. The disturbed profile of the X- 
ray emission seen i n the ROSAT obse rvation suggests 
an ongoing merger ([Tucker et al.lll998l ). Furthermore, 
it is a ssociated with a powerful radio halo ([Liang et al. I 
l200lf ), probably radiating via the synchrotron process, 
which suggests the presence of a population of ultra- 
relativistic particles. RX J0658-5557 was first observed 
by Chandra in 2000 Oct ober with 24.3 ks of usable data 
([Markevitch et afll2002[ ), and now has a total of over 500 
ks of Chandra exposure clearly revealing the bow shock 
in front of the "bullet" emerging from the merger. The 
collision is clearly supersonic, with a Mach number of 3 
deduced from the angle of the Mach cone. 

The XMM-Newton data for the c l uster have 
been analyzed previous ly in I Zhang et al. I ([2004) and 
iFinoguenov et~atl ([2005h . who find an average temper- 
ature of 13.6 ±0.7 keV, using data in the 2-12 keV range 
and fixing nn to the Galac tic value of 6.5 x 10 20 cm -2 
([Dickev fc Lockman I I1990D . In addition, the XMM- 
Newton data have been used in joint spectral fits with 
the Rossi X-Ray Timing Explorer (RXTE) data (in 
the context of t he se arch for hard X-ray emission) by 
IPetrosian et al.l ([2006) , where the hard X-ray flux would 
be by Compton scattering of the cosmic microwave back- 
ground by the same relativistic particles that are re- 
sponsible for the ra dio emission. In their analysis, 
IPetrosian et al.l ([2006) use MOS and pn data in the range 
1.0 - 10.0 keV, and, adopting the absorbing column to 
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be that measured for our Galaxy of 4.6 x 10 cm 
by iLiang et al. I (|200lf ), they infer an average temper- 
ature for the region within a 4' radius of 12.0 ± 0.5 
keV, in marginal agreement with that determined by 
iFinoguenov et al.1 (|2005l ). with the difference probably 
being due to the use of different bandpasses, calibration 
files, and the assumed absorbing column density. 

Similarly to Abell 1689, the spatial structure of the 
gravitational potential in this cluster has been exten- 
sively studied using both strong an d weak gravitationa l 
lensing. Weak-lensing analysis by I Clowe et al.l (|2004l ) 
and, m ore recently, j oint w eak- and strong- lensing anal- 
ysis bv iClowe etall (|2006l ) and lBradac et al.1 (|2006l ). re- 
veals a striking offset between the peaks in the gravi- 
tating material and X-ray-luminous matter. This sug- 
gests a scenario in which the subclusters have collided 
head-on and the gas is being slowed down by ram pres- 
sure, while the dark matter is able to pass through more 
or less without resistance. As such, this cluster offers 
one of the most compelling arguments for dark matter, 
which interacts only via gravitation, being distinct from 
the baryonic material, which is responsible for the X-ray 
luminosity. 

2.3. The Centaurus cluster 

The Centaurus cluster, also known as Abell 3526, is one 
of the most nearby clusters, and because of this, it en- 
ables good spatially resolved studies of cluster structure. 
It is a relatively cool c luster, at an average t emperature of 
kT = 3.68±0.06 keV (jFukazawa et al.|[l9 98). with a mod- 
est luminosity. One of its distinguishing characteristics is 
the detection of spatially resolved gradients of elemental 
abundances, even with data of modest angular resolution 
(see, e.g., Fukazawa et al. 1994, 1998). This, as well as a 
wide distribution of temperat ures of the emitting plasma, 
is clear in the Chandra data (Fa bian et al.ll2005f ), where 
the X-ray image shows a "swirl" in the central struc- 
ture, with a filament extending towards the north-east. 
This feature could have been shaped by a strong mag- 
ne tic field or a bulk flow wit hin the intracluster medium. 
In landers fe Fabian! (|2006T ). the XMM-Newton data are 
analyzed in combination with the Chandra data in or- 
der to derive maps of abundances for separate elements. 
They find that the element ratios are consistent with the 
solar ratios, with a metallicity of up to twice the solar 
value. A recent observation using the XIS detector (the 
X-ray Imaging Spectrometer) aboard the Suzaku satel- 
lite finds no evidence for bulk motion of the cluster gas 
and put an upper limit of |Ai;| < 1400 km s - 1 on any 
line-of-sight velocity difference (|Ota et al.ll2007t ). 

2.4. Observation details and data reduction 

The data were reduced using standard pipeline pro- 
cessing as of XMM-Newton SAS version 6.5 producing 
photon event lists. For the screening of soft proton flares, 
we create light curves in the 10-12 keV band for MOS 
and in the 12 - 14 keV band for pn in 100 s bins. We 
discarded the data when the total flux reached 3 a above 
the quiescent level (cf. iPratt fc Arnaudl l2002h . After 
this cut, we perform a similar second cut on the soft 
flux light curves in the 0.3-10 keV band for MOS and 
the 0.3-12 keV band for pn, binned by 10s. We use all 
event patterns (singles, doubles, triples and quadruples) 
for MOS and singles and doubles only for pn. We also 



TABLE 1 

Observation parameters 



Name / ObsID Detector Exposure (ks) Count Rate (s 1 ) 



A 1689 

0093030101 MOS 37 4.4 

pn 29 15 

RX J0658-5557 

0112980201 MOS 25 3.3 

pn 22 9.9 

Centaurus 

0046340101 MUS 46 18 

pn 42 60 



require XMMEA.EM for MOS and XMMEA.EP for pn, and 
also FLAG=0. The resulting effective exposure times and 
gross count rates in the 0.3 — 10 keV band for the three 
observations are shown in Table [TJ 

In the case of A 1689, we use t h e sam e data set as re- 
ported in lAndersson fc Madejskil (|2004l ) but reprocessed 
with the latest software and calibration data. For the 
"bullet" cluster, we extracted and analyzed data col- 
lected during the XMM-Newton pointing on 2000 Octo- 
ber 20 - 21; this i s the same XMM-Ne w ton observation as 
was report ed inlZhang et al. 1 (|2004fh IFinoguenov et ail 
(|2QQ5f ) and iPetrosian et al.1 (|2006f ) fsee above). Finally, 
the Centaurus observation was collect ed on 2002 January 
3 and is the same as the one used bv lSanders fc Fabianl 
(|2006h . 

Using the SAS command eexpmap, we create exposure 
maps with 1" x 1" bins in detector coordinates for all 
detectors. These account for bad pixels and columns in 
the data and correct for varying exposure over the CCDs. 

3. APPLICATION OF SPI TO XMM-NEWTON DATA 

3.1. Smoothed Particle Inference 

The method of Smoothed Particle Inference (SPI), 
which was constructed to model diffuse X-ray-emitting 
astrophysical sources, is described in PMA07. Here we 
only give a brief summary of the basic features of the 
method. 

In order to describe currently observed diffuse X-ray 
sources, a model with thousands of parameters is re- 
quired. Our choice of model is a set of spatially Gaussian 
X-ray emitters of an assigned spectral type. Each Gaus- 
sian "blob" is described by a spatial position, a Gaussian 
width, and a set of spectral parameters. 

In the Monte Carlo simulation, the flux of the astro- 
physical source at each energy and spatial position is con- 
verted to a prediction of the number of photons detected 
at a given detector position and energy. We calculate 
the probability of detection, D, given the instrument re- 
sponse R and the source model flux F: 



D(x,y,p) = j dEdO d0 R{x,y,p\0,^E) — — — . 

(1) 

Here (x, y) is the position on the detector, p is the ob- 
served pulseheight, (j)) are sky coordinates and E is 
the photon energy. This integral is calculated by sim- 
ulating photons sequentially while taking into account 
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mirror and detector characteristics (described in the next 
section, 3.2). 

The goodness of fit of the model is calculated from the 
likelihood function of the model parameters. The model 
data consist of a finite number of simulated photons, 
and we use a two-sample likelihood statistic to assess 
the goodness of fit. We explore the parameter space of 
the model with the Markov chain Monte Carlo (MCMC) 
method. The high dimensionality of the parameter space 
requires a method that is capable of exploring this space 
without being trapped in local minima. The Markov 
chain step is Gaussian, with a width that varies depend- 
ing on parameter history. We also find it necessary to 
make some parameters coupled. These are the same for 
all particles (global) and vary simultaneously. 

3.2. The XMM-Newton EPIC response function 

The ESA XMM-Newton satellite consists of three co- 
aligned X-ray teles copes and an optical/UV telescope 
(|Jansen et al.ll200H ). The telescopes focus X-rays onto 
two reflection grating spectrometers (RGS) and onto 
three CCD arrays for imaging spectroscopy: MOS1, 
MOS2 and pn. The imaging detectors are collectively 
referred to as the Europe an Photon Imaging Camera 
(EPIC, iTurner et al.lE_00lh . 
The full details of EPIC are covered in lEhle et aTl 
and its latest calibration is described in Kirsch 



2006). This section briefly describes the EPIC response 



function as calculated by us using the XMM-Newton 
Current Calibration Files (CCF) library 1 and explains 
how it is used in the Monte Carlo. A detailed descrip- 
tion of the response calcu l ation and interpolation is de- 
scribed in iPeterson etHI (|2004f l for the XMM-Newton 
RGS Monte Carlo. F or EPIC the detector response is 
different in structure, but the Monte Carlo process is the 
same. 

In summary, the response probability function of the 
EPIC cameras can be written as 



R(x,y,p\0,<i>,E)=A(E)v( 



,E)psf(x,y\6,^E)v R (E) x 

f(E)q(E)r(p\x,y,E)i(x,®] 

where A(E) is the mirror effective area, v(0, </>, E) is 
the mirror vignetting, psf(x, y\6, 0, E) is the energy- 
dependent point-spread function (PSF), vr(E) is the 
Reflection Grating Array (RGA) vignetting, f(E) is 
the filter transmission, q(E) is the quantum efficiency, 
r(p\x, y, E) is the pulse-height response, and i(x, y) is an 
exposure map correcting for bad pixels and differences in 
exposure time between different CCDs. 

This function gives the probability of detecting a 
photon of energy E originating at sky coordinates (0, (j)) 
as an event of pulse height p at detector coordinates 
(x,y). The details of these parts of the detector response 
are given below. 

Effective area and vignetting: The effective mir- 
ror collection area of the XMM mirrors is a function of 
energy and decreases with off- axis angle. This decrease 
is known as vignetting and is caused by shadowing from 
neighboring mirror shells. We obtain the on- axis effective 
area A(E) and the vignetting v(0, 0, E) for each telescope 



1 All c alibration data were obtained 
ftp: / /xmm. vilspa.esa.es/pub/ccf/constituents/ 



from 



from the XMM calibration files. 2 The response files are 
interpolated linearly and rebinned in order to optimize 
the performance of the Monte Carlo simulation. The bin 
size for A(E) is set to 50 eV, ranging from to 10.45 
keV. The v(6, 0, E) bins are 0.75 keV wide in energy and 
0.01° wide in off- axis angle. 

Point-spread function: We use a circular symmet- 
ric approximation for the PSF as it is described in the 
calibration files. For on- axis extended sources of < 10' 
in radial extent, this is a sufficient approximation for our 
purposes. Here we approximate the PSF by use of the 
encircled energy for different photon energies available 
from the CCF. 3 We use encircled energy samplings for 
photon energies from to 9 keV with a 1.5 keV spacing. 
For each energy band the encircled energy is sampled at 
1.37" intervals out to 9' from the center of the PSF. 

RGA vignetting: In the MOS cameras a fraction 
of the photons are intercepted by the RGA with some 
energy dependence. This loss of photons for MOS1 and 
MOS2 is tabulated in the calibration files. 4 In our Monte 
Carlo model, we use a sampling of this energy depen- 
dence of 0.5 keV. 

Filter transmission: The transmission of the optical 
blocking filters aboard XMM-Newton is modeled for the 
thin, medium and thick filter configurations as a function 
of energy. We rebin the transmission as found in the 
existing calibration files 5 into bins of 50 eV. 

Quantum efficiency: The ability of the detector 
CCDs to detect photons as a function of energy, or quan- 
tum efficiency (QE), is tabulated in the CCF response 
files. 6 We rebin the quantum efficiency for the Full Frame 
mode for MOS and both the Full Frame and Extended 
Full Frame modes for pn. We bin the QE in the range 
from to 12 keV, with 15 eV spacing. 

Pulse height redistribution function: EPIC pn 
response matrices are available from the XMM-Newton 
SAS Web site. 7 In the pn detector, the CCD response 
of the 12 CCDs varies with the distance from the line 
separating the two CCD rows. The pn response matrices 
are available for every 20 pixel rows from rows 0-20 (Y0, 
at the edge) to rows 181-200 (Y9, at the center). In the 
Monte Carlo simulation, we calculate the distance from 
the detector center line and use the correct pn response 
matrix accordingly. We only use matrices for the Full 
Frame and Extended Full Frame observing modes, and 
only for event patterns 0-4 (singles and doubles) for pn. 

EPIC MOS response matrices are dependent on the 
observing epoch and should be chosen according to the 
satellite revolution when the observation was taken. We 
use the SAS command rmf gen to generate MOS response 
matrices for 14 different epochs from revolution 101 to 

2 XRT1_XAREAEF_0008.CCF XRT2_XAREAEF_0009.CCF 
and XRT3_XAREAEF_0010.CCF. 

3 XRT1_XENCIREN_0003.CCF, XRT2_XENCIREN_0003.CCF 
and XRT3_XENCIREN_0003.CCF. 

4 RGS1_QUANTUMEF_0013.CCF and 
RGS2_QUANTUMEF_0014.CCF in FITS extension 
RGA.OBSCURATE. 

5 EMOS1_FILTERTRANSX_0012.CCF, 
EMOS2_FILTERTRANSX_0012.CCF and 
EPN_FILTERTRANSX_0014.CCF. 

6 EMOS1_QUANTUMEF_0016.CCF, 
EMOS2_QUANTUMEF_0016.CCF and 
EPN.Q UANTUMEFj0016.CCF. 

7 See ftp:/ /xmm. vilspa.esa.es/pub/ccf/constituents/extras/responses/ 
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revolution 1021. In the Monte Carlo we choose whichever 
epoch is closest to the observation. For MOS we use 
imaging mode matrices with all event patterns (0-12; sin- 
gles, doubles, triples and quadruples). Recently, it has 
been discovered that the MOS response is also dependent 
on distance from the detector axis (see the XMM-Newton 
EPIC Response and Background File Page, update 2005- 
12-15 8 ). In the XMM CCFs, the response is modeled 
in three different regions: a "patch," "patch wings," and 
outside "patch." Therefore, we also generate response 
matrices for all three regions for each epoch. In the MC 
simulation all three are read in for the epoch in question, 
and the correct one is chosen on the basis of the location 
of the detected photon. 

The pn response matrices are rebinned to an 800 x 800 
matrix with a constant 15 eV binsize from 0.05 to 12.05 
keV for both energy and pulseheight. MOS matrices have 
a 15 eV binsize and run from to 12 keV. These matri- 
ces are integrated over pulseheight to form cumulative 
distributions. 

3.2.1. Response Calculation 

In the Monte Carlo, photons are generated via proba- 
bility density functions normalized to unity. We calculate 
the cumulative distribution by integrating the probabil- 
ity functions and then draw a number from to 1 at 
random to choose a particular photon property. First, 
photons are chosen from a model function with a given 
set of input model parameters. The output variables are 
the photon energy and sky coordinates (E,9,(j)). Sec- 
ond, we predict the detector coordinates and pulseheight 
(p, x,y) by drawing photons using a response function 
R. In order to maintain the proper effective area and ex- 
posure normalization, photons are sometimes discarded 
according to the proper response functions (i.e. mirror 
effective area, filter transmission, vignetting, quantum 
efficiency, and exposure map). The response functions in 
equation ([2]) are in general not analytic and have to be 
stored in memory on grids. In order to limit the amount 
of used internal memory, we save the functions on coarse 
grids and interpolate linearly to get intermediate values. 

3.3. Application to the Data 
3.3.1. Data Binning 

In the framework of SPI, a three-dimensional adaptive 
binning technique is used to bin the photon event lists. 
The details of this are described in PMA07. The bins 
are rectangular in shape and are created so that each bin 
with 20 photons or more is split in two starting out from 
the entire dataspace. The choice of which dimension (x,y 
or pulse height) should be split is random, with a pos- 
sibility to split one dimension more often, on average, 
than another. Here we choose to split the pulseheight 
dimension on average 10 times more often than the x 
and y dimensions. This results in three-dimensional bins 
that are fine in pulseheight and coarse in x and y, suit- 
able for spatially resolved spectroscopy. The sizes of the 
bins are defined with respect to the full size of the data 
space in that dimension (20' for x and y and 9.7 keV for 
pulseheight). 

8 See http://xmm.vilspa.esa.es/external/xmm_sw_ca 1/calib/ 
epic_nles.shtml 



For RX J0658-55 this results in 14361 bins with 13.3 
photons bin -1 on average, and for the A1689 data we get 
31096 bins which equals 13.5 photons bin -1 on average. 
For the Centaurus data, the binning resulted in 197092 
bins with 13.7 photons bin -1 on average. 

For analysis, we use MOS data in the 0.3 - 10.0 keV 
range and pn data in the 1.1 - 10.0 keV range. The pn 
low-energy data are c ut off due to pn-MOS disagree ment 
at low energies (cf. lAndersson fc Madejskil l2004f ). In 
both analyses we restrict ourselves to a 20' x 20' spatial 
region centered on the respective cluster centers. After 
the model photons have been generated, they are binned 
on the same grid as the data photons before the two- 
sample likelihood is calculated. 

3.3.2. Spatially Resolved Spectral Model 

To model the cluster emission, we use a multi- 
component model consisting of spatially Gaussian 
smoothed particles, or "blobs," of cluster emission. Each 
of these is described by a spatial position, a Gaussian 
width, a single temperature, a set of elemental abun- 
dances, and an overall flux. Since each particle is de- 
scribed by a Gaussian, there will always be overlapping 
components. This means that the model everywhere de- 
scribes a multi temperature plasma. 

In this case we set the spectr al emission model 
to b e described by the MEKAL dMewe et all 119851 . 
119861 : lKaastralfl992l : iLiedahl et al.lll995f ) thermal plasma 
model with solar abundances absorbed by a WABS 
(M orrison fc McCammonl fl983l ) absorption model. The 
prior ranges for all parameters are set to be flat for a 
fixed range, except for the spatial Gaussian sigma which 
has a logarithmic prior distribution. The midpoint of 
the spectral parameter ranges is determined from sim- 
ple spectral analysis using the full cluster emission. The 
width of the range is chosen from the values that are ex- 
pected for that parameter in the cluster. It is, in general, 
an advantage to choose a parameter range that is wider 
than that expected. We choose to let both the equivalent 
hydrogen column and the redshift of the cluster plasma 
be variable in the analyses, as well as temperature and 
metallicity with respect to solar. The ranges used for the 
different clusters are shown in Table O 

In order to describe both the spatial and spectral prop- 
erties of the clusters adequately, we choose to use 700 
particles for the A1689 analysis and 600 for RX J0658- 
55 and Centaurus. A justification for this number of 
particles for data sets of similar complexity can be found 
in PMA07. However, we also do an analysis using only 
100 components in order to check the consistency of the 
broad spatially varying spectral properties of the clus- 
ters. 

3.3.3. Background Model 

The XMM-Newton instrumental background consists 
of three different main parts: particle background (soft 
protons), cosmic-ray-induced internal line emission, and 
electronic noise. In our background model, epicback, 
the particle background is approximated by a power-law 
spectrum with a variable spectral index and a separate 
normalization for each detector. This part of the model is 
not propagated through the mirror model, but is exposed 
directly onto the CCDs. In fact, this background compo- 
nent is scattered somewhat by the mirrors and does have 
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TABLE 2 
Allowed parameter ranges 



Parameter 


Min value 


Max value 


Global? 


A 1689 


n H (10 2U cm 2 ) 





3.5 


Y 


T (keV) 


5 


11 


N 


Z/Zq 


0.15 


0.35 


N 


z 


0.15 


0.21 


N 


R.A. & decl. a (arcmin) 


-10' 


+10' 


N 


In a (arcsec) 


0.5" 


5.5" 


N 


RX J0658-5557 


n H (10 2U cm^) 


2.7 


3.7 


Y 


T (keV) 


1 


19 


N 


Z/Zq 


0.12 


0.32 


N 


z 


0.28 


0.30 


N 


R.A. & decl. a (arcmin) 


-io' 


+10' 


N 


In cr (arcsec) 


0.25" 


4.25" 


N 


Centaurus 


n H (10 2U cm^) 





2.2 


Y 


T (keV) 


0.5 


9.5 


N 


Z/Zq 


0.1 


1.5 


N 


z 





0.015 


N 


R.A. & decl. a (arcmin) 


-10' 


+10' 


N 


a (arcmin) 


0.25" 


4.25" 


N 



a W.r.t. the nominal pointing of XMM. 

a radial dependence (Re ad et "aT1l2005f ). It decreases to 
about 80 % of its central value 10' from the pointing axis. 
This effect will be included in future papers. However, 
here, since we are dealing with bright clusters, we as- 
sume that a spatially flat modeling of this background is 
sufficient. 

We determine the best-fit parameters of our back- 
ground model using several EPIC observations with the 
filter wheel in the closed position 9 and "blank fields" 
with removed point sources, 10 as well as background files 
compiled by the XMM-N ewton Science Operations Cen- 
tre (SOC; fl^m^l2002h . We find the best-fit value of 
the power-law index to be approximately —0.22 (varying 
from -0.20 to -0.24 in the different observations), and we 
choose to fix it at this value in the remaining analysis. 

The internal line emission constitutes of fluorescent 
lines excited by high energy particles in various mate- 
rials of the detector. These lines include the Al K, Si K, 
Au M, Cr K, Mn K, Fe K, Ni K, Cu Ka, Cu K/3, Zn 
K, Au La and Au L/3 complexes. The lines are approx- 
imated with delta functions at the respective energies, 
and this is a good approximation, considering the lim- 
ited energy resolution of the CCDs. The emission is as- 
sumed to be uniform across the detectors, except for the 
Cu K emission, which is highly nonuniform, with a hole, 
devoid of emission, in the center of the detector plane. 
We approximate this hole with the intersection of a cir- 
cle with a 390" radius and a rectangle 630" wide. The 
relative normalization of these lines is determined indi- 
vidually for the MOS and pn detectors, using the filter 
wheel-closed observations. 

The electronic noise background includes bright pixels 

9 Observation IDs 0073740101, 0134521601, 0134522401, 
0134720401, 0136540501, 0136750301, 0150390101, 0150390301, 
0154150101, 

0160362501, 0160362601, 0160362801, 0160362901, and 
0165160501 

10 Observation IDs 0147511601 and 0037982001 



and columns, readout noise, etc., and is modeled as an 
exponential, F oc e~^ E / Ei \ with a variable value of Ei, 
where E is the energy assigned to the noise event. We 
find an appropriate value for E{ to be ~ 150 eV when 
using the MOS cutoff at 0.3 keV and the pn cutoff at 1.1 
keV. This value is a fixed parameter in the analysis. This 
noise background is assumed to be spatially uniform. 

The fraction of photons going to each of these model 
components (particle, lines, and noise) is variable in the 
analysis and all have a flat prior from to 1. Also, the 
fraction of photons for each model component going to 
each detector (MOS1, MOS2, and pn) is variable. The 
total normalization of the background model with re- 
spect to the cluster model, as well as the hard and soft 
X-ray background (XRB) components, is then set as a 
variable parameter also with a fiat prior from to 1. 

We model the soft Galactic X-ray background using 
a uniform emission component consisting of a MEKAL 
spectral model with WABS absorption. The plasma 
temperature is fixed at 0.16 keV, with a metal abun- 
dance of O.3Z at z = 0, and the absorption is fixed 
at nn = 1.5 x 10 20 cm -2 . Similarly, the hard XRB, 
presumably due to superposition of unresolved AGNs, 
is modeled using a power law with T = 1.47 and with 
absorption fixed at n# = 1.2 x 10 21 cm -2 . In general, 
it is customary to use zero absorption for the soft XRB 
and Galactic a bsorption for the hard comp onent in XRB 
analysis (e.g. iHickox fc Markevitchl (20 06) . We choose 
here to use the above values simply due to the fact that 
they give a better fit to the data in our analysis of the 
blank fields mentioned above, as well as source- free re- 
gions outside the clusters analyzed here. 

An example of the background model spectrum for two 
observations with the filter wheel in the closed position 
is shown in Figure [TJ All parameters of the background 
models are global. 

3.3.4. Markov Chain Model Sample 

In the Monte Carlo, photons are simulated according 
to the probability functions given by the model parame- 
ters, as described in Section UTTl The simulated photons 
are propagated through the detector model and binned 
on the grid determined by the three-dimensional photon 
density of the data, as described in Section 13.3.11 The 
two-sample likelihood Poisson statistic is calculated to 
provide a goodness of fit. In our analysis we use a model- 
to-data oversimulation factor of 10 for all data sets to 
reduce the model noise. In principle, it would be ideal 
to utilize a factor as large as possible, but we are limited 
by finite CPU speed and internal memory. In PMA07 
we compare different values of this factor and show that 
the results improve significantly when using a value of 
10 or greater. Parameters are iterated by Markov chain 
sampling, as described in the previous Sections. 

Figure [2] shows the evolution of the statistic (Poisson 
X 2 /dof.) — 1. Illustrated are (from left to right) values 
for the 700 blob run for A1689, the 600 blob run for 
RX J0658-55 and the 600 blob run for Centaurus. The 
statistic (Poisson x 2 /dof) approaches a value very close 
to 1 and stabilizes after ~ 2000 iterations. The value 
of the statistic at iteration 2000 for the three data sets 
(1.005, 1.013 and 1.016) is taken as an indicator that we 
have achieved an accurate fit. 

In the bottom panels of the same figure we show the 
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Fig. 1. — Typical background model including electronic noise, internal line emission, and particle background from a filter wheel-closed 
observation in pn (left) and for both MOS detectors (right). The flux is shown in units of counts per 25 eV, with the data shown as dashed 
lines and the model as solid lines. The ratio of the two is shown in the lower panels. The effective exposure time for the pn exposure is 28 
ks, whereas for each MOS exposure it is 23 ks. 



evolution of the (Poisson x 2 /dof)-l value for the control 
runs with 100 smoothed particles. In general, the runs 
with more blobs reach a good fit faster and approach a 
lower value of the fit statistic. On the basis of these plots 
we consider the chain to be stable after 2000 iterations, 
when the statistic becomes stationary and close to 1. We 
only use samples from iteration 2000 on to deduce cluster 
properties in the following sections. 

4. FINAL MODEL 

To visualize the output model sample we use the meth- 
ods described in PMA07. We stop the Markov chains 
after 4000 iterations and assume convergence after 2000 
iterations, as described above. We are left with 2000 
models, each consistent with the data. Each iteration 
takes approximately 30 minutes on a single Intel Pen- 
tium 4 2.0 GHz CPU which results in about 3 months of 
computing time per cluster. The models are filtered and 
marginalized in order to deduce the cluster properties 
that are discussed below. 

4.1. Abell 1689 

To confirm that we have an acceptable overall spec- 
tral fit, we plot the model spectrum as inferred from the 
model sample versus the data and the ratio of the two in 
Figure 03 This plot shows an overall accurate spectral fit 
with the exception of the instrumental Cu K line, which 
appears to be detected at a slightly higher energy than 
expected. This mismatch is most likely due to gain vari- 
ations in the XMM-Newton CCDs. The accuracy of the 
fit is seen by looking at the statistic in Figure El 

To infer the luminosity distribution of the cluster, we 
generate the luminosity map for each model in the sam- 
ple and take the median in each spatial pixel in order 
to make the median luminosity map shown in Figure [H 
(right). The bolometric luminosity, L^, is displayed in 
units of 10 44 erg s _1 per 2" x 2" . The map of raw counts 
(counts per 2" x 2") from the screened data file is shown 
on the left for comparison. There is good agreement be- 
tween the reconstruction and the raw data; specifically, 
the reconstructed profile is sharper than the data due to 
PSF deconvolution and the fact that obvious chip gaps 



in the data are compensated by taking the exposure into 
account. 

We have selected three interesting regions (num- 
bered in Figure IH), motiv ated by our analysis in 
lAndersson fc Madejskil ([2004) to study in more detail 
the distribution of plasma temperatures in those regions. 
First, for a cluster that has a quite regular surface bright- 
ness distribution, similar to clusters with well-established 
"cooling cores," the temperature of the core for Abell 
1689 is quite high: ~ 7.5 keV. Second, we have identified 
a region of hotter plasma, ~ 9 keV, to the north of the 
cluster core, indicating possible shock heating of the gas 
in that region. Third, we have chosen to study in detail a 
region south of the core exhibiting a temperature almost 
as low as that of the core itself. 

Next, we form a median temperature map by first 
weighting the temperature of each model particle by its 
luminosity, creating a temperature map for each model 
in the 2000 model MCMC sample. This sample of maps 
is then averaged by taking the median of the distribu- 
tion of temperatures for each spatial pixel, as is done in 
the creation of luminosity maps. Instead of taking the 
luminosity- weighted average, in another attempt to visu- 
alize the distribution of temperatures in the cluster, we 
bin the luminosity in bins of temperature in each spatial 
pixel, creating a three dimensional differential luminos- 
ity data cube. We do this for all sample models, and for 
each of them, we calculate the mode of the resulting dis- 
tribution of each spatial pixel. Taking the median over 
all model samples of this mode accurately describes the 
dominant temperature at any given spatial position. We 
note that this will not give the same value of the temper- 
ature as one would get from traditional spectral analysis. 
The temperature measurement is biased due to the de- 
generacies present in spectral fitting when allowing for 
700 separate temperature phases. However, we believe 
that it is the method of highest contrast when separat- 
ing distinct temperature components of a galaxy cluster. 
The median distribution mode and the median emission- 
weighted temperature are shown in Figure [5] (top left 
and top right) . To display an estimate of the error on the 
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Fig. 2. — Poisson x 2 per degree of freedom for (top, from left to right), the Abell 1689, RX J0658-5557, and Centaurus standard runs 
(700,600 and 600 blobs respectively), and (bottom, from left to right) the Abell 1689, RX J0658-5557, and Centaurus 100 blob runs . 



emission- weighted temperature we calculate the 1 a vari- 
ation on these quantities over the model sample. These 
maps are shown in the same figure (bottom left and bot- 
tom right). The uncertainty of the emission- weighted 
temperature is below 0.5 keV throughout most of the 
shown region. For all of the calculations above we have 
used the model run with 700 particles. 

In order to study the regions selected above in more 
detail, we have extracted the differential luminosity dis- 
tribution in these regions and binned them into 20 bins 
over the 5-11 keV allowed range. We have chosen to use 
the model run with 100 particles for this in order not 
to over-complicate the problem. While in this case the 
luminosity and temperature maps do not show the same 
level of detail, the distributions of spectral parameters 
become narrower. 

The differential luminosity, plotted as the emission 
measure, j n#n e dV, per keV, along with the median 
of the distribution, is shown for these regions in Figure 
[6] From these data alone we cannot detect any devia- 
tion from isothermality in these regions, mostly due to 
the inherent degeneracies associated with the attempt to 
fit high-temperature emission with a multi-phase model. 
Since most of the emission is from the bremsstrahlung 
continuum there are many combinations of gas phases 
(temperatures) that provide equally good fits to the data. 

We have confirmed the existence of a region of hot 
plasma north of the c luster core as we reported in 
lAndersson fc Madejskil (|2Q04f ). This hotter plasma ap- 
pears to extend in an arc approximately halfway around 
the northern part. The technique verifies both the pres- 
ence of the colder emission to the south and the lower 
than ambient (by ~ 1 keV) temperature of the cluster 
core. Our analysis reveals an apparent trail to the south, 
possibly a remnant tracing the path of the cluster core, 
and the shock-heated region to the north. However, this 
could also be due to projection effects from filaments ex- 
tending in the direction of the line of sight. 

In our earlier findings we infer intra-cluster gas mo- 



tions from a shift in the position of the Fe K line com- 
plex corresponding to Az ~ 0.01. Here we show a map 
of cluster redshift (Figure [7J left) that is similar to the 
temperature maps shown above. The map was created 
by calculating the dominant redshift mode at each spa- 
tial pixel and taking the median over the samples. The 
superposed contours are iso-photon count from the raw 
data, and the grid is shown for comparison with the re- 
gions us^dJnjAnd^ (|2004l ) . The results 
from lAndersson fc Madejskil (|2004l )~are shown in Figure 
El (right). The map clearly confirms the Az ^0.01 east- 
west gradient found previously by us, corresponding to 
a velocity of ~ 3000 km s _1 . We also assess the signifi- 
cance of this by analyzing the redshift-mode variation in 
the model sample. The la model variance of the redshift 
mode in the regions of interest ranges from 0.016 to 0.019, 
and therefore our result is only ~ la significant. Putting 
an upper threshold on the redshift discrepancy, we find 
that Az < 0.045. However, we note that it is not clear 
that the above is the most appropriate estimate of the 
true uncertainty, since inherent properties of the method 
could increase the variance; possibly, another means of 
error estimation could give a more accurate result. Also, 
this particular model was not designed specifically to de- 
duce the redshift-structure of the cluster gas. Devising 
the run only to measure this effect will likely prove more 
successful. 

4.2. RX J0658-55 

Similarly to the case of Abell 1689, we plot the model 
spectrum as inferred from a subsample of the model sam- 
ple versus the data and the ratio of the two in Figure [3l 
This plot shows an acceptable fit, again with the excep- 
tion of the Cu K line complex. The evolution of the 
Poisson x 2 as a function of Markov chain iteration is 
shown in Figure [2j As described previously, we generate 
a median luminosity map from the model and compare it 
with a counts map smoothed by a 4" Gaussian (Figure [8] 
right and left). It can be seen that the "bullet" and the 
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Fig. 3. — Spectral comparison of the data (dashed line) to the model (solid line) for the Abell 1689 700 blob analysis (top left), the RX 
J0658-55 600 blob analysis (top right), and the Centaurus 600 blob analysis (bottom) for the full 20' x 20' field used in the analysis. Flux 
is in units of photon counts per 25 eV bin, with the ratio of data to model shown in the lower panels. The model spectra are produced by 
averaging the spectrum from every 100th iteration from 2000 to 4000 and are renormalized to match the data counts. The spectra include 
all backgrounds and all three EPIC detectors. 
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0.004 0.006 



Fig. 4. — Results of the analysis for Abell 1689, showing the central 5' x 5' region. Left: Raw data smoothed by a A" kernel Gaussian 
(counts per 2" x 2" pixel). Right: Luminosity reconstruction using the median of all samples at each spatial point (Lfo o ^/10 44 erg s _1 per 
2" x2"). 

main cluster features become somewhat sharper by the PSF deconvolution effect from the forward fitting that is 
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Fig. 5. — Results of the analysis for Abell 1689, showing the central 5 X x 5 ; region. Top left, temperature map (in units of keV) constructed 
using the mode of the distribution for all samples at each spatial point; top right, temperature map (in units of keV) constructed using the 
median of the distribution for all samples at each spatial point; bottom left, temperature uncertainty map (in units of keV) showing the 1 
a variation of the temperature mode in the model sample; bottom right, temperature uncertainty map (in units of keV) showing the 1 a 
variation of the emission-weighted temperature in the model sample. 
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Fig. 6. — Iteration averaged distribution of temperatures for A1689 regions 1 (core; left), 2 (south; middle) and 3 (north, right). 



inherent in the method. 

We form and plot the luminosity weighted temperature 
map as well as a map based on the dominant temperature 
mode, in Figure [9] (top right and top left). In conjunc- 
tion, the associated uncertainty maps are shown (bottom 
right and bottom left) . Both maps clearly show the cold 
core remnant of the "bullet" and the hot shocked gas in 
front of it. Note that the mode and weighted tempera- 
ture maps display different properties of the temperature 



structure and are not comparable. This of course is due 
to the approach used to model the plasma in our method. 
Instead of using just a single temperature, as is custom- 
ary, in every spatial region, we use a number of phases 
that over the model sample form a nearly continuous dis- 
tribution. We discuss this further in the paragraph dis- 
cussing isothermal simulations below. The mode is often 
biased towards lower temperatures for plasmas with tem- 
peratures of 9 keV and lower, due to the fact that within 



11 




Fig. 7. — Left: Redshift map of Abell 1689 obtained by taking the median of the dominant redshift mode over the model sample. The 
gray scale from black to whit e corresponds to z = 0.175 — 0.185. Superposed contours from the raw event list data a nd gr id regions from 
Andersson & Madeiski (2004) are shown for comparison. Right: Results from the analysis in Andersson &; Madeiski (2004) 



the bandpass of EPIC, at higher T, the spectra are more 
similar to each other and therefore are more difficult to 
distinguish. 

We do not see a clear decrease in temperature in front 
of the shock, as is expected, since the gas should be undis- 
turbed here. It is possible that emission from the post- 
shock gas is smeared by the XMM-Newton PSF and com- 
pletely dominates the pre-shock emission, which in turn 
would be very faint in comparison. The cold gas of the 
bullet shows an apparent tail stretching south-east from 
the bullet center. By looking at the temperature map 
alone, one would conclude that this tail might reveal the 
movement history of the bullet in the merger. However, 
the luminosity map and, much more clearly, the 500 ks 
Chandra exposure (M. Markevitch et al. in preparation), 
show a symmetrical Mach cone directed westward, indi- 
cating that this is the direction of motion. The weak- 
lensing analysis (|Bradac et al.ll2006f ) also shows the west- 
ern dark matter halo just west of the bullet. It still can- 
not be ruled out that the bullet core entered to the north- 
east of the main cluster (see the "entry hole" devoid of 
emission in this region), proceeding south-west through 
the main cluster and eventually slingshotting around to 
the northwest, thus creating the tail visible in the tem- 
perature maps. The dark- matter dominated regions and 
the gas-dominated ones, of course, do not need to have 
similar merger paths. This scenario is also supported by 
the apparently previously shock-heated region of gas lo- 
cated directly south of the cluster. This part seems to 
have almost as hot a gas as the shock in front of the 
bullet. We select this region (No. 6) along with the six 
other regions shown in Figure to do a more detailed 
analysis. 

In Figure [10] we show the detailed distributions of tem- 
peratures, in numerical order, from the selected regions 
above. In region 1 we have extracted the central emis- 
sion from the cluster core. The distribution clearly shows 
the signatures of a kT = 7 keV plasma, with possible 
contamination from higher temperatures (seen as peaks 
around 11 and 14 keV), most likely from projection ef- 
fects. Regions 2-7 show the distributions of various re- 



gions in the cluster. 

In order to determine the deviation from isothermal- 
ity of the plasma in the selected regions we generated 
isothermal data sets of a cluster resembling RX J0658- 
55 with the same number of photons and the same spatial 
structure. The isothermal models were reconstructed us- 
ing the exact same method as in the reconstruction of 
the real data. The distribution of temperatures from the 
isothermal reconstructions are shown for 4, 7, 10, and 
15 keV plasmas in Figure [Til None of the distributions 
in Figure [10] can conclusively be distinguished from an 
isothermal plasma. 

Finally we have created luminosity maps of RX J0658- 
55 using emission components in separate bands of tem- 
perature. Figure [12] shows the cluster luminosity in the 
kT = 1 — 7, 7 — 13 and 13 — 19 keV bands. The contours 
are iso-luminosity contours from the kT = 1 — 7 keV 
map. These maps show an interesting property hinted 
at by the temperature maps: the bullet appears to move 
further north with higher temperature. It is possible 
that this feature indicates an increased compression to 
the north-west resulting in a higher temperature. This is 
supported by the Chandra image which shows a shorter 
distance between the bullet and shock front to the north- 
west than to the southwest. 

The three-phase map also clearly shows the exten- 
sion of hotter gas to the south that is causing the high- 
temperature region in the temperature maps (region 6). 

4.3. The Centaurus Cluster 

The model spectrum of the Centaurus cluster as in- 
ferred from the model sample is shown with the data 
spectrum and the ratio of the two in Figure [3] The model 
spectrum can be seen to fit the data well, and this is a 
consequence of the fact that we have a low value of the 
overall Poisson x 2 , as can be seen in Figure [2] However, 
there are some residuals at > 5 keV energies due to some 
inadequately fitted spectral lines. Since most of the data 
are in the low-energy range, these data tend to drive the 
fit, leading to some unwanted excesses and deficits at 
higher energies. These, however, do not have a major 
impact on the statistic. 
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Fig. 9. — Results of the analysis for RX J0658-55, showing the central 5 X x5 x region. Top left, temperature map (in units of keV) constructed 
using the mode of the distribution for all samples at each spatial point; top right, temperature map (in units of keV) constructed using the 
median of the distribution for all samples at each spatial point; bottom left, temperature uncertainty map (in units of keV) showing the 1 
a variation of the temperature mode in the model sample; bottom right, temperature uncertainty map (in units of keV) showing the 1 a 
variation of the emission-weighted temperature in the model sample. 

We form both luminosity and temperature maps, anal- [13] the raw count map of the XMM-Newton data (left) 
ogous to our procedures in previous sections. In Figure is shown as well as the luminosity reconstruction (right). 
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Fig. 10. — Iteration-averaged distribution of temperatures for RX J0658-55 in regions 1-7, also showing the median. 




Fig. 11. — Iteration-averaged distribution of temperatures for the reconstruction of isothermal simulations of RX J0658-5557 due to 4, 
7, 10, and 15 keV plasmas. 



Even though the gaps from dead pixel rows are taken 
into account via the exposure map, some artifacts can 
be seen. The model aligns itself with the chip gap where 



a filament extends north-east of the cluster core. Figure 
[Ml shows a temperature mode map (top left), as well as 
a median temperature map (top right), with the relevant 
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Fig. 12. — Luminosity maps (L^/10 44 erg s 
the ranges are 1 - 7, 7 - 13 and 13 - 19 keV. 



per 2" x 2") of RX J0658-55 in three different temperature bands. From left to right, 



uncertainty maps shown below. The dominant tempera- 
ture mode can be seen to be elevated with respect to the 
median in a hot region north-east in the direction of the 
filament. 

To investigate the temperature in some of the more 
interesting regions in detail, we calculated the temper- 
ature distributions in four regions and compared them. 
These plots are shown in Figure [15l The order of the 
regions corresponds to the numbers in Figure HU Se- 
lected regions correspond to the cluster core (No. 1), 
the extended filament (No. 2), the ambient temperature 
directly to the west of the core (No. 3), and the anoma- 
lously hot region to the north-east (No. 4). The core 
shows signs of nonisothermality; the distribution includes 
a narrow peak at 0.5 keV, as well as a bump around 1.5 
keV, very probably due to projection. In the other plots 
it is harder to distinguish the distributions from isother- 
mal. However in region 4 there is a hint of a hotter 
~ 8keV phase. 

For the Centaurus cluster we have focused on a new 
approach in analyzing cluster structure. In Figure [TBI we 
have produced median maps for the luminosity in bands 
of different temperatu res. This is similar to the method 
of IFabian etHI (|2006h . who fit a six-phase plasma model 
to finely spatially binned data of the Perseus cluster in 
order to make maps of gas mass in six different tempera- 
ture bands. The main differences with our modeling are 
that here we use a smooth cluster model and allow for a 
very large number of variable temperature phases. From 
left to right, Figure [16] shows the cluster luminosity in the 
- 1, 1 - 2, 2 - 4, 4 - 6, 6 - 8, and 8 - 10 keV temperature 
bands. 

This subdivision reveals some striking features: the 
cluster core is dominant in the 0-1 keV band, and ap- 
parently it has moved in from the northeast, leaving a 
tail of colder gas. The 1-2 keV temperature band map 
shows a bi-polar nature of this emission around the core. 
With higher temperatures, the emission becomes more 
and more offset, and the 8-10 keV map shows a con- 
centration of superheated gas in an isolated region to the 
north-east, possibly the remnant of a merger. This was 
no ted in the analys i s of th e Chandra data for this cluster 
by I Crawford et all ([2005). Interestingly, this feature is 
aligned with a filament extending from the cluster core. 



Possibly these two irregular phenomena are related. 

Finally, we explore the metallicity structure of Cen- 
taurus by creating a median metal abundance map, 
analogous to our construction of median temperature 
maps. The metallicity map is shown in Figure [T71 and 
confirms previous findin gs by IFabian et al.l ([2P05) and 
ISanders fc Fabianl (|2006[ ) that the metallicity increases 
by about a factor of 2 (from 0.5 to 1.0 solar) toward s 
the center. However, in contrast to IFabian et al.1 (|2005l ). 
we do not find metal abundance levels as high as twice 
the Solar value and we also do not see a sharp decrease 
towards in the very center. A minor dip can how- 
ever, be distinguished. It is most likely that both the 
absence of very high metallicity and the sharp central fea- 
ture are due to PSF smearing effects. We also confirm a 
tight correlation between temperature and metallicity in 
the gas, with lower temperature plasma generally having 
a higher metallicity. Calculating the Pearson product- 
moment correlation coefficient for the median values of 
temperature and metal abundance per 2" x 2" bin, as 
shown in Figures [TH and [T71 weighted by the luminosity 
in that bin gives a coefficient of —0.87, hence confirming 
a strong correlation. We do not find any significant dis- 
crepancy in the spatial distribution of th e redshift of the 
gas, in agreement with lOta et al.l (|2007l ). 

5. DISCUSSION 

This paper describes an application of the Markov 
chain Monte Carlo technique developed by us for the 
analysis of observations with the XMM-Newton imaging 
instruments. We demonstrate the flexibility and power 
of this technique - employing smoothed particle infer- 
ence - via studies of three very different clusters, Abell 
1689, RX J0658-55 and Centaurus, especially regarding 
the ability to determine the spatial distribution of tem- 
perature of the radiating plasma, which is difficult via 
more traditional techniques. 

We found evidence for cluster merger activity in all 
these systems, but in each case, the signature was quite 
distinct. The bullet of RX J0658-55, the remnant of a 
merger in Centaurus and the asymmetry of temperature 
in A1689 may roughly correspond to the early, middle, 
and late stages of cluster merging. In all cases the core 
of the cluster seemed relatively unaffected. Further sys- 
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Fig. 13. — Results of the analysis for the Centaurus cluster, showing the central 6 f x 6 X region. Left: Raw data smoothed by a A" 
kernel Gaussian (counts per 2" x 2" pixel). Right: Luminosity reconstruction using the median of all samples at each spatial point 
(L 6oi /10 44 erg s" 1 per 2" x 2"). 
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Fig. 14. — Results of the analysis for the Centaurus cluster, showing the central 6 f x 6 ; region. Top left, temperature map (in units of 
keV) constructed using the mode of the distribution for all samples at each spatial point; top right, temperature map (in units of keV) 
constructed using the median of the distribution for all samples at each spatial point; bottom left, temperature uncertainty map (in units 
of keV) showing the 1 a variation of the temperature mode in the model sample; bottom right, temperature uncertainty map (in units of 
keV) showing the 1 a variation of the emission weighted temperature in the model sample. 



tematic studies in temperature structure may add to our ties. 

understanding of the effect of mergers on cluster proper- The most important difference of our technique com- 
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Fig. 15. — Iteration-averaged distribution of temperatures for Centaurus in regions 1-4 also showing the median. 




Fig. 16. Luminosity maps (L bo j/10 44 erg s _1 per 2" x 2") of the Centaurus cluster in the - 1, 1 - 2, 2 - 4, 4 - 6, 6 - 8 and 8 - 10 keV 
temperature bands. The superposed contours represent full band luminosity. 



pared to conventional modeling is that we use overlap- 
ping emission components. This facilitates the use of 
a large number of phases describing the X-ray emis- 
sion at each spatial coordinate. This scenario is much 
more physical than when one adapts the assumption that 
any given point in the projected cluster image can be 
described by a single-temperature phase. However, it 



also poses the challenge of constraining a distribution of 
phases when the spectra of high-temperature plasmas in 
particular are very similar. Comparison with other anal- 
yses requires the projection of a distribution into a single 
value. We have chosen in this paper to show the mode 
and the median of that distribution in our temperature 
maps along with the associated 1 a variations (see Fig- 
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Fig. 17. — Metallicity map of the central 6' x 6 ; region of Centaurus showing the metallicity with respect to solar, and with contours 
representing luminosity. 



uresEJEl and[T4j). 

The difference in the assumptions in modeling pre- 
cludes any direct comparison with other measurements, 
and it is more than likely that the results will be different. 
The median of our distributions in particular will depend 
on the choice of the allowed range of temperatures in the 
model. This is especially true for the choice of the upper 
boundary in high-temperature sources. The selection of 
the boundaries must be based on prior knowledge of the 
source, such as measurements in the hard X-ray band. 
This introduces a systematic uncertainty that has to be 
taken into account. 

Another advantage of this technique is the way we treat 
the background. Since the background is modelled in- 
stead of subtracted, difficult aspects of background sub- 
traction, such as weighing the different vignetting of dif- 
ferent components, are avoided. The background model 
has been calibrated against a number of observations 
in which the filter wheel was closed, as well as against 
numerous "blank sky" observations (see Section f3.3.3|) . 
Even though we allow the normalizations of different 
background components to be optimized by the MCMC, 
there is still a small systematic uncertainty associated 
with residual calibration uncertainties. 

In this work we have opted to simultaneously fit the 
spatial distributions of the temperature, metallicity, and 
redshift of the X-ray-emitting, plasma illustrating the 
broad power of the method. This is not necessarily the 
optimal approach to answer a specific question about 
the cluster, such as resolving a redshift difference across 
A1689. This may explain the larger statistical uncer- 
tainty obtained in this work compared to o ur conven- 
tional analysis (|Andersson fc Madejskill2Q0l ). For such 
a problem the method could be designed specifically for 
this purpose. Since temperature and metallicity varia- 
tions are small across this cluster these could be global 
parameters in the modeling significantly simplifying the 
overall problem. For such tasks that are highly depen- 
dent on detailed spectral features, the number of cluster 
particles should also be kept to a minimum to avoid the 
risk of smearing the spectral feature in question. This is, 



however, a compromise, since a lower number of particles 
would limit the spatial resolution of the spectral feature 
that one wants to resolve. 

In summary, the statistical uncertainties in this work, 
which are inherent in problems with large numbers of 
degrees of freedom, may seem larger than those that are 
achieved in traditional analyses, but the uncertainties 
can always be reduced by simplifying the problem. This 
work attempts to constrain all cluster characteristics into 
a single model. There will be other applications in which 
a specific model can be tested to answer a specific prob- 
lem. 

6. FUTURE WORK 

We expect that this method will prove to be power- 
ful in the determination of cluster gas density and en- 
tropy. Since the model is multi parametric and analyti- 
cal, based on the Markov chain posterior, the smoothed 
particles can be manipulated in order to construct a 
three-dimensional cluster model. In principle, the z- 
coordinate can be chosen for each particle on the basis 
of some assumption of spherical symmetry. This can be 
done separately in ranges of particle plasma tempera- 
ture, assuming that particles of similar temperatures ex- 
hibit the same structure as in the two-dimensional case. 
This will accurately produce the temperature gradient 
in the three-dimensional case while preserving the two- 
dimensional observed structure. This is likely the most 
accurate method to determine the three-dimensional 
spectrally resolved structure of galaxy clusters in the X- 
ray band. 

Since the likelihood calculation does not have to be 
limited to a single instrument, but can potentially in- 
clude multiple X-ray instruments or even data from other 
kinds of measurements (such as gravitational lensing, 
Sunyaev-Zeldovich data, or optical velocity dispersions), 
the method is quite general. This method is also well 
suited to analyze other complex, spatially resolved ob- 
jects where the observed X-ray emission might consist of 
superposed separable components with varying spectral 
parameters, such as in supernova remnants. 
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